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ABSTRACT 

Starburst regions with multiple powerful winds of young massive stars and supernova 
remnants (SNRs) are favorable sites for high-energy cosmic ray (CR) acceleration. A 
supernova (SN) shock colliding with a fast wind from a compact cluster of young stars 
allows the acceleration of protons to energies well above the standard limits of diffusive 
shock acceleration in an isolated SN. The proton spectrum in such a wind-SN pevatron 
accelerator is hard with a large flux in the high-energy-end of the spectrum producing 
copious 7 -rays and neutrinos in inelastic nuclear collisions. We argue that SN shocks in 
the Westerlund 1 (Wdl) cluster in the Milky Way may accelerate protons to > 40 PeV. 
Once accelerated, these CRs will diffuse into surrounding dense clouds and produce 
neutrinos with fluxes sufficient to explain a fraction of the events detected by IceCube 
from the inner Galaxy. 

Key words: neutrinos — acceleration of particles — ISM: cosmic rays — ISM: super¬ 
nova remnants — magnetohydrodynamics (MHD) — shock waves 


1 INTRODUCTION 


Gamma-ray observations with both space- and ground-based 
telescopes have shown that supernova remnants are the main 
sources of CRs at least up to 100 TeV. Galactic sources of PeV 


CRs, however, are still to be identified (see e.g. Blandford 
|et al.|[^014[ |Amato|[2014[ ) although it has been argued that 
type Ilb supernovae, a subclass which comprise about 3% of 
the observed core collapse SNe, may be able to produce CRs 


with energies beyond 100 PeV (e.g., Ptuskin et al. 20101. 


Cosmic neutrino observations are well suited for identifying 


galactic pevatrons (see e.g. Halzen & Hooper 2002 Aharonian 
2004 Becker 2008[ Anchordoqui et al. 20141. The IceCube 


South Pole Observatory (IceCube) has detected 37 neutrino 
events above the expected atmospheric neutrino background 
in a 988-day sample ( Aartsen et al.|[2014a l, with three PeV 
neutrinos, I.O4I+J44 PeV, 1.Idling PeV and 2.004t262 PeV, 
being the most energetic neutrino events in history. While 
significant spatial or time clustering of the events has not 


* E-mail: byk@astro.ioffe.ru 
f E-mail: don_ellison@ncsu.edu 
t E-mail: peter.gladilin@gmail.com 
§ E-mail: osm2004@mail.ru 


yet been reported, a possible association of some events with 
galactic center sources was proposed (e.g., Razzaque|2013 1. 


The ANTARES neutrino telescope | Adrian-Martinez 
et al. 12014 1, using six years of data collected near the galac¬ 


tic centre, reported 90% confidence level upper limits on the 
muon neutrino flux to be between 3.5 and 5.1 xlO~® GeV 
cm“^ s“^, depending on the exact location of the source. 
They excluded a single point source as the origin of 7 neu¬ 
trinos observed by IceCube in the vicinity of galactic centre. 
However, an extended source of a few degrees is not excluded. 

Since the most likely high-energy neutrino producing 
mechanisms are the inelastic p—nuclei and p — 7 collisions 
of protons, where the reaction kinematics result in the en¬ 
ergy of the neutrinos to be ~ 0.05 that of the protons, 
the energy of the parent protons should exceed 4 • 10^® eV 
to explain the IceCube observations. The 7-rays produced 
in these reactions have ~ 0.1 of the proton energy (e.g., 
Halzen|[2013 1. Furthermore, the proton accelerators must be 


very efficient to produce the high-energy neutrino flux of 
i/Fi, « 10“® GeV cm“^ s“^ sr“^ per flavor in the 0.1-1 PeV 
range detected by IceCube. 

Neutrinos from photo-meson p — 7 interactions in com¬ 
pact particle accelerators, like the cores of active galactic nu- 
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clei (AGN) (e.g., Stecker 2013 \ and 7 -ray bursts (GRBs) (e.g., 


Waxman fc Bahcall|1997 Meszaros fc Waxman|200l| , along 


with other models (e.g., Ahlers & Murase 2014| Fox et al. 


2013 

[Kashiyama et af. [2013 He et af. 2013 

Murase et al. 

2013 

[Liu et al.|2014| Kashiyama & Meszaros 

2014| Tambor- 


to explain the origin of these first-ever IceGube neutrinos. 


It has been suggested by Neronov et al. (20141 that Ice 


Gube neutrinos and Fermi/LAT 7 -rays are both produced in 
interactions of GRs which have a hard spectrum with a power- 
law index harder than 2.4 and a cut-off above ~ 10 PeV. This 
assumes the GRs are interacting with the interstellar medium 
in the Norma arm and/or in the galactic bar. 

The role of isolated SNRs in producing PeV GRs is un¬ 
certain. The 7 -ray spectra from most shell-type galactic SNRs 
observed by Gherenkov telescopes show a spectral cutoff well 


below PeV energies (see e.g. Acero et al. 20151. While fu¬ 


ture Cherenkov Telescope Array observations may be needed 
to confirm these results (see e.g. Aharonian|2013 (, a popula¬ 
tion of pevatrons is needed to explain the observed spectrum 
of galactic GRs at and above the spectral knee region. Su¬ 


pernovae with transrelativistic shocks (Budnik et al. 2008 


Chakraborti et al. 20111 and type Ilb supernovae (Ptuskin 
et aI.||20T0| were proposed to accelerate GRs above PeV en¬ 


ergies but the statistics of these potential sources remain to 
be established. Here, we present a model of a galactic peva- 
tron which produces hard CR spectra with a high efficiency of 
conversion of SN kinetic energy into the highest energy GRs. 
We further argue that this PeV CR source has the proper¬ 
ties required to explain a number of the IceGube neutrinos 
detected from the direction of the inner Galaxy. 

Core-collapse supernovae in associations of massive OB 
stars are certain to produce a fraction of galactic GRs, as 
demonstrated by isotopic measurements by the Advanced 


Composition Explorer (see e.g., Binns et al. 20071. Active 
star-forming regions may comprise extended associations of 


massive young stars like Cyg OB2 (see e.g., Wright et al. 


20141, as well as compact dense clusters of young massive 


stars like Westerlund 1 (Wdl). These extended and compact 
cluster types are distinct in both spatial and temporal scales 
(e.g., Gieles fc Portegies Zwart]|2011 1 . Both type of clusters 
of massive stars are expected to be efficient PeV CR acceler¬ 
ators in starburst and normal galaxies (see Bykov|[2014 for 
a review). 

The compact massive cluster Wdl, with an estimated 
age of 3.5—5 Myr, contains more than 50 post main sequence 
stars (e.g., Clark et al. 2005| , including at least 24 Wolf-Rayet 
(WR) stars of both flavours (representing about 8 % of the 
observed galactic population of WR stars) in about a parsec- 
scale core (e.g., Clark et al. [2005 Crowther et al. 20061. 
These massive stars have strong individual winds with an 
estimated total cluster kinetic power exceeding 10 ^® erg s“^. 
In the small, dense cluster core these individual winds should 
combine and drive a fast cluster-scale wind by the mecha¬ 
nisms studied by |ChevaIier fc Cle^ (|1985[) and [Stevens "fcl 


Hartwell (20031. The magnetar CXOU J1647-45, discovered 


by Muno et al. 


(2006b I using high resolution Chandra X-ray 


was likely produced about 10,000 years ago by a supernova 
with a progenitor star of mass ^ 40 Mq (Muno et al.|200^ 


Mereghetti|2008 1 and remains the only direct evidence of su¬ 


pernova activity in Wdl. 

In our pevatron model, a SN blast wave collides with the 
termination shock of a strong wind generated by the collec¬ 
tive action of many massive stars in a compact cluster. Both 
shocks are assumed to propagate in a homogeneous upstream 
plasma. We show that proton energies well above a PeV may 
be produced with a hard spectrum where the CR spectrum 
at PeV energies is most sensitive to the shock speeds and am- 
plihcation of magnetic turbulence associated with CR driven 
instabilities. The model provides a high efficiency of conver¬ 
sion of SN shock-cluster wind ram pressure into PeV GRs 
which is needed to explain at least some of the IceGube neu¬ 
trino events. We don’t consider a SN interacting with its own 
wind since individual stellar winds will be unimportant in 
the compact cluster environment where dozens of massive 
stars are located within a few parsecs and the wind cavity is 
smoothed out on this scale. 

We show that diffusive shock acceleration (DSA) in sys¬ 
tems with colliding shock flows (CSFs) can provide maximum 
particle energies, CR fluxes, and energy conversion efficien¬ 
cies well above those produced in an isolated SNR shock of 


the same velocity (Bykov et al. 20131. While CSFs are ex¬ 
pected to occur in colliding stellar winds, the most powerful 
events should happen when a supernova shock impacts the 
extended fast wind of a nearby young massive star cluster. 

The important features of shock acceleration in CSFs 
are: (i) the production of a piece-wise power-law particle dis¬ 
tribution with a very hard spectrum of confined particles at 
the high-energy end just before a break, and (ii) an increase 
in the maximum energy of the accelerated particles, and the 
acceleration efficiency, compared to that obtained with DSA 
at an isolated SNR shock of the same speed. These two prop¬ 
erties imply that a substantial fraction of the flow ram pres¬ 
sure is converted into relativistic particle pressure. The high- 
energy GRs, therefore, must modify the dynamics of the CSF 
system. To model the spectra of accelerated particles in CS¬ 
Fs, a nonlinear, time-dependent model was constructed in 
[Bykov et al.| ( [2013[ ) . The maximum energy and absolute flux¬ 
es of the accelerated particles, both inside the CSFs, and 
those escaping the acceleration site, depend on the shock ve¬ 
locities, the number densities, and the magnetic fields in the 
flows. We show below that sub-PeV and PeV neutrinos from 
colliding shocks in galactic and extragalactic compact clus¬ 
ters of young stars with reasonable parameters, such as those 
expected in the Wdl compact cluster, can explain a fraction 
of the IceGube neutrino events. 


2 PARTICLE ACCELERATION IN COLLIDING 
SHOCK FLOWS 

In order to model the proton acceleration in a compact stellar 
cluster, we used the nonlinear, time-dependent model of CS- 


observations, has been associated with Wdl. This magnetar 


Fs presented in Bykov et al. (20131. The main modification 


is that here we allow for different parameters for the wind 
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Figure 1. Colliding flow geometry where the stellar wind (left) 
was approximated from the analytic model of ( |Wilkin|1996^ , while 
the SNR shock was assumed to be spherical. The star, or cluster 
center, approaches the SN explosion center with speed Vq- The 
spectrum of the Fermi accelerated CRs was derived taking into 
account effects of the flow velocity projection on the curved shock 
surfaces. 


termination shock and the SNR blast wave. We have also in¬ 
troduced an approximation to account for the shape of the 
shock (as illustrated in Fig. using an analytic expression 
by Wilkin (19961. Other approximations, such as the Bohm 


diffusion in the shock vicinity, and a parametrization of the 
magnetic held amplification (MFA) due to CR-driven insta¬ 


bilities, are the same as in Bykov et al. (20131. The reader 
is referred to that paper for full details. Our results for neu¬ 
trino and 7 -ray production in Wdl, given in Section use 
our nonlinear, time-dependent model. However, to describe 
the general characteristics of CSFs, we shall start in this sec¬ 
tion with simple analytic estimates from a linear model in a 
plain-parallel case (e.g., Bykov et al.||20lT K 

Considering a SNR expanding in a compact OB- 
association; at some expansion phase the distance L between 
the SNR blast wave and a stellar wind shock is less than the 
mean free path of the highest energy CRs in the SNR shock 
precursor. At this point, the CR distribution function around 
the two shocks, indicated by i = 1,2, can be approximated 
as 


p{x,p,t) = Ap ®exp (-^ |a;|) X 
xH (p-po) H {t- ts.cc) , 


where the CR acceleration time is 

V 

tacc — 

PO 


[ 3 / Di dp 

J {Ui + U 2 ) \Ui U2 J p 


( 1 ) 


( 2 ) 


ft is important to note that, besides only applying to 
high-energy particles with mean free paths larger than L, 
these equations are qualitatively different from those of test- 
particle DSA in an isolated shock: the spectrum fi below 
the exponential break is harder and the acceleration time is 
shorter. 

Our model assumes a high level of CR-driven magnetic 


instabilities and Bohm diffusion for CRs in the close vicinity 
of the shocks (see, e.g., Bell||2004 Schure et al.||20f2 Bykov 


et al.|2014 1. However, once high-energy particles obtain mean 


free paths on the order or larger than the distance between 
the shocks, L, scattering will become much weaker. A specific 
feature of the simulation is that the highest energy CRs, with 
Pmax ^ P ^ P*, propagate with little scattering. Despite the 
weak scattering between the shocks, their momenta are still 
nearly isotropic since they scatter for long periods in the ex¬ 
tended regions downstream from the shocks. Here, p* is the 
momentum such that the proton mean free path A(p*) > L. 
Since the particle distributions are nearly isotropic even for 
high-energy particles with A(p) > L, the kinetic equation re¬ 
duces to the so-called telegrapher equation (e.g., Earl||1974 1 
and this allows a smooth transition between the diffusive and 
the scatter-free propagation regimes. 

The time-dependent nature of our simulation means that 
protons escape the accelerator at different stages of the sys¬ 
tem evolution (i.e., as the SNR blast wave approaches the 
stellar wind termination shock) producing pions and neutri¬ 
nos with varying hardness and maximum energy and these 
effects are included in our results. Of course, the complex evo¬ 
lution of the source and some unknown details of the mass 
distribution in the outer ISM region (e.g., the presence of 
dense shells or clouds) will also influence the results. How¬ 
ever, we believe the general properties of our simulation are 
robust. 

Assuming Bohm diffusion with D\{p) = D 2 (p) = D(p) — 
cRg{p)/3, due to CR-driven, amplified magnetic instabilities 
in the CR accelerator, one obtains 


cRgjp) 


( 3 ) 


where ui = Us is the SNR shock velocity, U 2 = u™ is the 
stellar wind speed, and Rg{p) is the momentum dependent 
proton gyroradius. Then, using the scaling B « ^/Ajvrjbpui 
for the amplified magnetic field, the acceleration time can be 
estimated as face ~ 2 • 10^° epev {Vbn)~^'^ u~f u~\ s. In the 
above expressions, p = mpU, rjb is the acceleration efficiency, 
i.e., the fraction of ram kinetic energy in the plasma flows 
converted into accelerated particles, the energy is in PeV, 
and the speeds are in units of 10^ km s“^. 

While the ejecta speed in the free expansion SNR 
phase can have a wide range of values, we use a mean 
ejecta speed of ~ /Mq km s“^ and take 

the mean duration of the free expansion phase to be ~ 
200{M/Mq yr. In the specific case of a young 

SN shock propagating through the winds of massive stars in 
cluster Wdl, where n ~ 0.6cm“^ (see 
, and assuming rjb ~ 0.1, we find face ~ 400 yr 
accelerated to epev ~ 40, when Us 3 ~ 10 and 
Um 3 ~ 3. The high SN shock velocity Us 3 ~ 10 is expected in 
the free expansion stage if the fast ejecta mass is about one 
solar mass. 

Therefore, for SNRs with ages less than Isnr ~ 400 yr, 
one can get face < tsNR for 40 PeV protons with standard 
parameters. In this case, at tsNR ~ 400 yr, the SNR radius is 
~ 3 — 4 pc. Furthermore, the hard spectrum expected from 


the compact 
et al.|[2Cl()6b 
for a proton 


Muno 
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CSFs puts most of the energy into the highest energy pro¬ 
tons and one can estimate the power in the highest energy 
neutrinos produced in inelastic p—p-collisions by the decay 
of charged pions (tt^ —^ eVe Vfi i'll) as 


(l038 cmO (s.OOOkms-i) (iP^) ® 

where is the fraction of energy in the inelastic p — p- 
collisions which is deposited in the high-energy neutrinos, S 
is the cross section of the colliding flows, and Tc is the con¬ 
finement time for protons in the emission region where the 
target density is Js 4 times the ambient density n due to 
the shock compression, Vg. 

In Eq. , we take the inelastic p—p collision cross section 
to be ~ 70 mb above a proton energy of 10 PeV. Then the pro¬ 
ton cooling time, nt^y, can be estimated as ntpp « 2.5 • 10^"^ 
s assuming that two inelastic collisions are needed to con¬ 
vert most of the proton energy into secondaries. The fraction 
of the proton energy deposited into neutrinos of all types, 
fv ~ 0.15, was derived using both the analytical parameter- 
izations for the inelastic p—p-collisions presented by |Keln^ 
|et al.| p006[ ); [i^fexhiu et al.| ( |2014[ ) and the GEANT4 pack¬ 
age simulations. At the energies of photons and neutrinos 
considered in Figs. m and [5} the neutrino fluxes calculated 
with the [Kelner et ( 2006| ~^ or [Kafexhiu et al.| ( |2014[ ) cross 
sections differ by less than 20 %. 

We note that the observed CR energy density and the 
galactic SN statistics require that rjp > 10% over the lifetimes 
of isolated SNRs to power galactic CRs. However, the instan¬ 
taneous efficiency may be much larger during early stages of 


the SNR evolution (Blasi 20131. In the case of a SNR colliding 


with a wind, we expect the efficiency to be even higher than 
in a young isolated SNR and assume it can reach Pp > 0.5 for 


the CSF stage lasting for a few hundred years ( Bykov et al. 
|2013[ ). Furthermore, since CSFs produce very hard spectra 
when face <C tesc, i.e., N(j) oc 7 “^, most of the CR energy 
lies in the high-energy tail just below the upper break which 
occurs when tacc is greater than the escape time, teac. The 
essential properties of CSF acceleration in young stellar clus¬ 
ters are high overall acceleration efficiency, spectra harder 
than N(j) oc 7 ”^, and maximum proton energies > 40 PeV 
achieved in a few hundred years. 


2.1 Colliding flows and 3D geometry effects 

While an exact treatment of non-linear CSF with a significant 
backreaction from accelerated CRs is unfeasible in 3D, we 
have introduced an approximation to account for effects from 
3D geometry in our plane-parallel model. Instead of taking 
the cluster wind termination shock and the SNR blast wave 
as plane, we account for aspects of the curved shock surfaces 
at positions parameterized by the angles 0 or </> in Fig.[2 We 
still assume planar shocks for the DSA calculation but with 
varying projected velocities at distances d(0) > L 12 along the 
shock surfaces away from the symmetry axis. Here, L 12 is 



Figure 2. The small dashed circles (not to scale) indicate the inner 
CR acceleration region in the massive young cluster WdlThe inner 
CR acceleration region has a radius of 3-4 pc (shown in a Chandra 
image) and the blow-up of the inner region shows the position of 
the magnetar CXOU J1647-45 found by |Muno et al.||2006b^ with 
a high angular resolution Chandra observation of WdlThe outer 
region is 30-40 pc centered around Wdl and indicated by the large 
dashed circle. This is overlaid on a H.E.S.S. map of TeV emission 
adapted from [Abramowski et al!] | |2012a'l l. The angular resolution 
of IceCube is larger than the outer circle. In Fig. we show the 
larger neutrino emitting ISM volume of radius ~ 140 pc around 
WdlThe volume is filled over ~ 10^ yr with CRs accelerated during 
a short CR acceleration phase ~ 400 years right after the supernova 
explosion in Wdl which produced the magnetar. 


the time-dependent minimnm distance between the colliding 
shocks. For each position determined by d(0) for 0 < 0 < 
90° we calculate the non-linear particle distribution using 
the projected speeds 14, cos 0 and Knr cos >1? as parallel flow 
speeds. That is, we assume the wind and the SNR shocks 
are locally plane with converging parallel flows set by the 
projected speeds. 

For this approximation, we nse the bow shock wind mod¬ 
el of Wilkin (1996 1 for R(0) and 14,, assume the SNR shock 
is spherical, and restrict our calculations to 3 arc surfaces: 
0 — 30°, 30 — 60°, 60 — 90° which produce accelerated particle 
densities with the weights 0.84, 0.15 and 0.01, respectively. 
Then the weighted CR distribution function is used to cal¬ 
culate the 7 -ray and neutrino emissivities both in the accel¬ 
erator (see (3.11 and in the surrounding ISM from the CRs 
escaped the accelerator (see ^3.2|. 
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Figure 3. Model predictions of 7 -rays (solid curves) and neutrinos 
(dashed curves) from p—p-interactions calculated in a CSF source 
of age 400 yr. The dotted curve is the inverse Compton emission 
from primary and secondary electrons accelerated directly in this 
source. The extreme upward curvature in the neutrino spectrum 
above ~ lOTeV reflects the transition from CR acceleration in 
the single SNR shock for low-energy particles to the more efficient 
acceleration for high-energy particles as they scatter back and forth 
between the SNR shock and the cluster wind. The data points 
for the H.E.S.S. source, and the five Ice Cube events explained 
in Fig. are presented to illustrate how they compare to our 
model predictions when the source is ~ 400 yr old and point-like, 
i.e., Wdl about 10“^ years ago. The simulated 7 -ray and neutrino 
emission at the present time from CRs that escaped the accelerator 
in Wdl ~ lO'^ yr ago are summarized in Figs. |^and[^ 

3 NEUTRINOS AND GAMMA RAYS FROM 
THE WESTERLUND 1 CLUSTER 

3.1 Emission from Trapped CRs 

To explain the IceCube neutrinos, we combine the nonlinear, 
time-dependent model of particle acceleration in CSFs with 
a propagation model applied to the Wdl compact cluster 
(Fig-i- The time-dependent simulations provide the evolv¬ 
ing energy spectra of CR protons and electrons as the SNR 
shock approaches the strong wind of a nearby early-type star. 
For relativistic electrons/positrons we account for the energy 
losses due to synchrotron and inverse Compton (IC) radia¬ 
tion. 

The particle acceleration is combined with a propaga¬ 
tion model relevant to Wdl. As described in Section i high- 
energy particles will obtain a hard spectrum when their mean 
free path is large enough so they scatter back and forth be¬ 
tween the two converging shocks. Lower energy particles how¬ 
ever, will be confined to, and accelerated by a single shock 
and obtain a softer spectrum. This phenomenon is illustrated 
in Fig. 3 of |Bykov et al.| ( |2013[ ). 

It is important to note that even though the CSF ac¬ 
celeration is efficient and the CR population modifies the 
structure of the plasma flows, the spectrum of high-energy 
CRs scattering between the two shocks remains hard until 


they gain enough energy to escape from the system. The 
shock modification can cause face to increase, and the energy 
where the exponential turnover in Eq. 0 starts to dominate 
drop, but below the turnover, the spectrum remains close to 
N{'y) oc 7 “^. 

Thus, the CSF system will have two spectral regimes for 
trapped CRs: a low-energy region (< 1 TeV) from particles 
accelerated in a single SNR shock (produced both before and 
after the start of the two-shock acceleration period), and a 
high-energy, hard spectrum region (> 1 TeV) from particles 
accelerated in the converging flows. The transition between 
these two regimes of acceleration occurs when the energy of 
a particle, Et, is large enough so it can easily travel between 
the two shocks. 

This transition, occurring at ~ 20 TeV, is seen as a bend 
in the neutrino spectrum shown in Fig. This figure shows 
neutrinos (solid curve) and y-rays (dashed curve) from p—p- 
decay from CR protons that are still trapped near the SNR 
and stellar wind, along with IC from trapped CR electrons. 
This is the emission expected ~ 400 yr after the SN explosion 
in a region of radius ~ 3 —4 pc. Note that the y-ray fluxes 
of the sources in Fig. are presented in erg cm“^ s“^ sr“^ 
in order to be compared with the observed diffuse neutrino 
fluxes, while the fluxes of the sources in Fig. are measured 
in erg cm“^ s“^ as usual. 

In Fig. we have overlaid a schematic of this inner CSF 
region on a map of Wdl. At later times, these “trapped” CRs 
will escape the accelerator and diffuse beyond the inner CSF 
region into a much larger outer ISM region where they may 
encounter dense clouds producing neutrinos and y-rays for 
an extended period of time, as we discuss in Section [3.2[ For 
clarity, the regions are not drawn to scale in Fig. 

To derive the gamma-emissivity of PeV CRs we account¬ 
ed for the Breit-Wheeler effect of pair production by energet¬ 
ic photons interacting with the interstellar radiation field as 


well as the extragalactic light background (see e.g. Aharoni- 
|an|2004j |Dwek &: Krennrich||2013l . This interaction leads to 
a significant suppression of the y-ray flux at the high-energy 
end of the spectrum for distant (^ 10 kpc) sources. 

For the emission from Wdl, with an estimated distance 
of ~ 4 kpc, this effect leads to a relatively small suppression 
of the y-ray flux at 1 — 10 PeV, as can be seen in Figs. 
and|^ In addition to the y-rays produced by pion decay, we 
have also included y-rays produced by the IC radiation of 
the secondary pairs produced in situ by the same p — p- 
interactions (see the IC curve in Fig. |^. In this calculation, 
we accounted for the pair energy losses in a mean magnetic 
field of magnitude 10 gG in the extended c loud of number 
density 25 cm“^ (see e.g. Strong et al.|20f4|. 


3.2 Emission from Escaping CRs 


The initial acceleration stage lasts a few hundred years after 
the SN explosion producing high-energy CRs that escape the 
CSF system and diffuse through the ambient ISM. However, 
the estimated age of the supernova which produced the mag- 


netar CXOU J1647-45 in Wdl is ~ 10 yr (Muno et al. 


2006a 


Mereghetti|2008 1. If the CSF acceleration occurred ~ 10 yrs 
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ago, these CRs will have produced pious over a much longer 
time span as the TeV-PeV particles diffuse away, fill a region 
of about 140 pc radius, and interact with the ambient ISM 
(the outer ISM region indicated in Fig. [^is about 30 pc). 

In Figs. 1^ and we show simulated 7 -ray and neutrino 
spectra produced by the CRs that escaped from the acceler¬ 
ator and diffused into the surrounding cloudy medium over 
10'* yr. At this point the acceleration responsible for the emis¬ 
sion in Fig. has ceased long ago. 

The diffusion model used to propagate the CRs has three 
regions. Within the accelerator of 3-4 pc radius we assume 
a Bohm diffusion coefficient Db = 3 x lO^^Fipev cm^ s“*. 
Outside of the dense cloud region we use an ISM value Db = 
3 X 10^®i5pp^ cm^ s“*, which is consistent with the standard 


models of CR propagation in the Galaxy (see e.g. Strong 


et al. 20071. Between these two regions we assume a transition 
coefficient Dtran = Ils(R/3pc)^, where R is the distance from 
the core. We assume a mean cloud density ~ 25cm“^, cloud 
radius ~ 30 pc, and have set, at 1 PeV, Db = Diqm at R = 
30 pc. A cloud of radius 30 pc with mean density 25 cm“^ 
would have a mass ~ 10® Mq. Outside the cloud we assumed 
an ISM density of ~ lcm“®. We note that the densities we 
assume for the clouds match available 7 -ray observations of 
Wdl ( [Abramowski et al.||2012^ |Ohm et al.|2013| ). 

The IceCube data points in Fig. show the neutrino 
energy flux consistent with the position of Wdl considering 
the angular resolution of the instrument. The CSF model 
applied to Wdl with reasonable parameters can explain a 
subset of the observed Ice Cube neutrinos. 



Figure 4. Gamma-ray emission from inelastic p— p-interactions 
in the CSF source at ~ 10* yr after the SN explosion when CR 
protons produced in the short-lived accelerator have propagated 
into a nearby cloud of ~ 30 pc size. The magnetic field amplified 
by the CR-driven instabilities in the vicinity of the fast shock in 
the CSF accelerator were parameterized as 0.8 mC (c), 0.9 mC (b), 
and 1 mC (a), all below 10% of the ram pressure. The 1C curve is 
inverse Compton emission from the secondary electrons produced 
by the inelastic p—p-interactions in the cloud. Only the 7 -rays from 
the H.E.S.S. field of view are included. The gas number density of 
the nearby cloud is 25 cm“®, except for the light-weight solid curve 
where it is 30cm“® with B = 0.9 mC. 


4 NEUTRINO PEV AND SUB-PEV EVENTS 


The median angular error values for the IceCube neutrino 
events are given in the Supplementary Material for |Aart-| 
sen et al. (2014aI. Using these median angular errors, and 
assuming Gaussian statistics (see Aartsen et al. 2014b| ), we 
have produced a sky map with 2 —cr contours (corresponding 
to about 86 % confidence in 2D Gaussian statistics) for the 
neutrino events in the vicinity of Wdl (see Fig. |^. As seen 
in this map, a source of neutrinos with a radius of ~ 140 pc 
around Wdl (black circle) can be associated with five neu¬ 
trino events, including two PeV events. These are 2, 14 (PeV 
event “Bert”), 15, 25 and 35 (PeV event “Big Bird”). 

Based on this sky map, we compared our model neutrino 
spectra from the Wdl source with the fluxes in five IceCube 
energy bins corresponding to the five neutrinos within 2 —ct 
of the source. This is shown in the Fig. A more precise 
comparison will require both more sophisticated models of 
the cloud distribution within a few hundred parsecs of Wdl 
and a more accurate determination of the event positions. 


5 DISCUSSION 

Explaining the origin of the recently detected PeV neutrinos 
by IceCube is a fundamental challenge for models of particle 
acceleration. The observations imply a source that can pro¬ 
duce substantial fluxes of protons with energies considerably 


higher than those expected from isolated SNRs. Since isolated 
SNRs remain the most likely source of CRs below a few PeV, 
the neutrino source must produce high-energy protons with¬ 
out conflicting with the observed properties of galactic CRs. 
The underlying protons producing the neutrinos in the CSF 
model described above have a hard spectrum, and are few in 
number, avoiding any conflict with low-energy CR population 
measurements. The CSFs may contribute to the high-energy 
end of the CR population produced by isolated SNRs and 
superbubbles (see e.g. Binns et al.||2007 1 and, we believe the 
strong plasma flows in compact clusters of young stars, such 
as Wdl, contain the energy and specific properties needed to 
explain a significant fraction of the IceCube neutrinos. 

Compact clusters contain massive stars with strong 
winds and recent SN activity. It is inevitable that occasions 
will occur when a SN blast wave collides with the termination 
shock from the strong wind of a nearby massive star, or with 
an extended cluster wind from several massive stars. We have 
developed a model of the Fermi acceleration expected from 
such colliding shock flows and, using realistic parameters, ob¬ 
tained simultaneous fits to the H.E.S.S. 7 -ray observations 
(Fig.|4f and the fraction of IceCube neutrinos expected from 
Wdl (Fig.[§. 
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Figure 5. Neutrino emission from an extended 140 pc radius) 
source ~ yr after the SN explosion when CRs produced in the 
short-lived accelerator have propagated into the surrounding ma¬ 
terial. The amplified magnetic fields are 0.8 mG (c), 0.9 mG (b), 
and 1 mG (a). We note that the 7 -rays in Figs. and the neu¬ 
trinos here originate from different volumes: only the 7 -rays from 
the H.E.S.S. field of view are shown in Fig. while the neutrinos 
are from a larger region of radius 140 pc. The neutrino data points 
(Icj energy flux error bars) are a subset from all 37 IceGube events 
| |Aartsen et al.|2014a[ |). These five events are within 2-a contours 
from Wdl based on Fig. Two PeV events (14, “Bert”) and 35 
(“Big Bird”), as well as three sub-PeV (2, 15, 25) events, are in¬ 
cluded in the subset. Note that the sub-PeV event 25 has a very 
large position uncertainty in Fig. [^and have to be considered with 
some care. 


5.1 Multi-wavelength signatures of CSF scenario 


A unique property of our CSF model is that the acceleration, 
while producing hard proton spectra to multi-PeV energies 
with high efficiencies, only lasts a small fraction of the SNR 
lifetime, just the time when the SNR is colliding with the 
nearby stellar wind. Because of the hard spectra of acceler¬ 
ated CRs most of the energy is in the highest energy regime. 
In Fig. 1^ we show the predictions for high-energy photon and 
neutrino emission of the source at the end of the accelera¬ 
tion stage (which was supposedly ~ years ago). The syn¬ 
chrotron radio emission then is estimated to be > 10 Jy at 

2.2 GHz, a value well above the current level. This is consis¬ 
tent with our scenario where the brief CR acceleration stage 
ended ~ 10 "^ years ago. 

Indeed, the total radio flux from Wdl as measured with 
the Australia Telescope Compact Array (ATCA) interferom¬ 
eter by Dougherty et al. (20101 is 422, 461, 523, and 669 mjy 


at 8 . 6 , 4.8, 2.2, and 1.4 GHz, respectively, and after subtract¬ 
ing the radio emission from stellar sources, they derived dif¬ 
fuse emission fluxes of 307, 351, and 426 mJy at 8 . 6 , 4.8, and 

2.2 GHz. Colliding winds in massive binary star systems were 


90 " 



Figure 6. Map showing 2-cr contours for a subset of IceCube 
events associated with the inner Galaxy as determined from 2D 
Gaussian statistics and the median angular errors and positions 
given by ( [Aartsen et al.|2014b[ |. Two PeV events (14, “Bert”) and 
35 (“Big Bird”) as well as three sub-PeV (2, 15, 25) events are 
within 2-cr from Wdl. 


proposed by [Eichler fc Usov| ( |1993| ) to accelerate relativistic 
particles and produce non-thermal radio and GeV regime 7 - 
ray emission. Some of the stellar radio sources detected in 
Wdl with the ATCA by [Dougherty et ah] ( |2010| exhibited 
composite spectra of both non-thermal and thermal emission 
potentially indicating particle acceleration in colliding wind 
binaries. 

In contrast to the binary wind model, CRs in CSFs are 
generated in the violent environment of a SN blast wave col¬ 
liding with a cluster wind. While the acceleration stage is 
brief, CRs will escape the source and interact with the near¬ 
by ISM clouds for long periods (> 10^yr). During this time 
pions will be produced in inelastic p—p-interactions resulting 
in relativistic secondary electrons and positrons from tt* de¬ 
cay and photons from tt® decay. The hard CR spectra from 
CSFs will result in prominent peaks in the spectral ener¬ 
gy distribution produced by synchrotron, IC, and pion-decay 
emission, as shown in Fig. For an extended cloud of size 
~ 30 pc and number density ~ 25 cm“^, the peaks correspond 
to keV, TeV, and PeV energy bands. 

For Fig. we assumed the cloud magnetic field to be 
B = 10 pG consistent with Zeeman splitting measurements 
of a number of molecular clouds compiled recently by jStrongj 


et al. (20141. Since the penetration of CR nuclei into the 
cloud may be reduced at low energies (e.g., jCesarsky &; Volkj 


1978 Protheroe et al. 20081, we show in Fig. [(^ spectra for 


the case with no proton penetration energy threshold (upper 
dashed curve) and for proton threshold energies F* = 10 
and 20 GeV. The radio fluxes are sensitive to F* but the 
X-ray synchrotron fluxes corresponding to the peak of the 
synchrotron SFD are not. 
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Figure 7. The spectral energy distribution of the synchrotron 
(dashed lines), inverse Compton (dot-dashed line) emission from 
secondary electrons and positrons, as well as photons produced by 
pion-decay (solid line) from the inelastic p —p-interactions in the 
nearby clouds which are H.E.S.S. sources shown in Fig. The 
CRs were accelerated at Wdl and diffused into the clouds. The 
cloud gas number density is 25 cm~^ with magnetic field B = 
10 pG. The 7 -rays detected by H.E.S.S. are indicated in the figure. 
The upper dashed synchrotron curve is the result with no energy 
threshold for CR protons to penetrade into the cloud. The lower 
synchrotron curves correspond to threshold values 10 GeV and 
20 GeV respectively. The threshold values are low enough to not 
influence the pion-decay emission. 


The cloud associated with the H.E.S.S. source in this 
model is a diffuse, low-surface-brightness (< 0.1 pJy arcsec“^ 
at 1.4 GHz), flat-spectrum, synchrotron source of polarized 
radio emission with bright spots of brightness ~ 1 mjy 


ments of the magnetic field in dense clumps. The future 
Square Kilometre Array, with a sensitivity of ~ 1 pjy/beam, 
may allow detection of radio emission from clouds irradiated 
by CRs (see e.g. Strong et al.|2014 1. 

The synchrotron peak in the CSF model is at keV X-ray 
energies and a source with a half-degree extension and to¬ 
tal flux of ~ erg cm“^ s“^ may be detectable with the 

future eROSITA (extended ROentgen Survey with an Imag¬ 
ing Telescope Array) instrument on the Spectrum-Roentgen- 
Gamma (see e.g. CappellutI] 20111 and ASTRO-H (see e.g. 
Aharonian et al. |2014| ) satellites . The search for synchrotron 


X-ray emission from the cold clouds located near these pow¬ 
erful CR sources may be conducted with the next generation 
of X-ray sky surveys. 

In Fig. we show simulated spectra of 7 -ray emission 
from CRs that escaped the accelerator and diffused into the 
surrounding cloudy medium over 10"* yr. The models are com¬ 
pared to H.E.S.S. data for the source associated with Wdl 


( [Abramowski et al.|2012a Ohm et al.|2013 1 and we have re¬ 
duced the 7 -ray flux to correspond to the H.E.S.S. field of 
view at Wdl. At 10^ yr, the acceleration responsible for the 
emission has long ceased and the emission comes only from 
CRs accelerated in the source and propagating through the 
ISM . 

We note that an analysis of 4.5 yr of Fermi-LAT data by 


Ohm et al. (20131 found extended emission offset from Wdl 


by about 1 degree. This study concluded that acceleration 
of electrons in a pulsar wind nebula could provide a natural 
explanation of the observed GeV emission. However, |Ohm| 
et al.| ||2013 1 found that the pulsar wind nebula could not 


explain the TeV emission observed by H.E.S.S.. As seen in 
Fig. a the CSF model can satisfactorily explain the TeV 7 - 
ray emission. 

There is an apparent excess of neutrino events, including 
two of the three PeV neutrinos in the IceCube map present¬ 
ed in ( [Aartsen et al.|[2014a (, within a radian from Wdl. In 
Fig. a we show this set in a map with the positions of the 
events and 2-a contours as determined from 2D Gaussian 
statistics and the median angular errors of the IceCube tele¬ 
scope. Westerlund 1 is indicated with a contour (black circle) 
corresponding to a 140 pc radius region - a few degrees - 
where neutrinos are produced in our CSF model as escaping 
CRs diffuse out from the compact accelerator for ~ 10'* yr. 
The neutrino energy flux corresponding to the solid curve in 
Fig. ais ~ 3.7 X 10“® GeV cm“^ s“* This is well below the 
90% confidence level upper limits imposed by ANTARES ob¬ 
servations ( Adrian-Martfnez et al.|2014 1 for a source of width 
> 0.5 degrees at the Wdl declination. 


5.2 CSFs in the starburst galaxies 


Another issue concerns how hard spectrum PeV neutrino 
sources contribute to starburst galaxy radiation. |Loeb fc] 
Waxman (2006 1 suggested that CR interactions in starburst 


galaxies may efficiently produce high energy neutrinos and 
contribute cumulatively into the neutrino background. The 
population of star forming galaxies with AGNs and the star- 
burst galaxies peaked at a redshift z > 1 , with a wide tail of 
the distribution revealed by Herschel ( Gruppioni et al.|2013 1 
up to z ~ 4.5. These objects most likely contribute to both 
the isotropic diffuse 7 -ray background measured by Fermi- 


LAT (Ackermann et al. 20151 between 100 MeV and 820 


GeV, and the diffuse flux of high-energy neutrinos measured 
by IceCube. 

Assuming that a CR spectral index for all the starburst¬ 
like galaxies is 2 . 1 - 2.2 at the high-energy part of the spec¬ 
trum, [Tamborra et al(] ( |2014[ ) were able to provide a rea¬ 
sonable fit to both the Fermi and IceCube data. Larger in¬ 
dices failed to explain the observed diffuse neutrino flux. 
That CR spectra harder than in normal galaxies like the 
Milky Way were required, may reflect a different population 
of CR sources and/or different CR propagation in the star- 
burst galaxies. Superclusters of young massive stars are likely 
much more abundant in starburst galaxies compared to the 
Milky Way since mergers and interactions of galaxies result in 
abundant supercluster formation (see e.g. [Conti et al.|2012[). 
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The hard spectra and high efficiencies from CSFs make them 
an attractive way to produce CRs well beyond PeV in star- 
burst galaxies where the collective contribution from many 
CSFs sources might extend the CR knee to higher energies 
compared to the Milky Way. 

High resolution radio observations of the star forming 
galaxies M82, Arp 220, NGC 253, M31, M33 and others pro¬ 
vide information on the magnetic field structure and lepton- 
ic CRs (Tabatabaei et ciL||2013 [Adebahr et al. 2013 Persic 


!& Rephaeli||2014[ ), and '^ay telescopes have observed some 

starburst galaxies up to ~ 10 TeV (see Acero fc et al.||2009 


Abramowski et al.|2012bl|Lacki et al.|201H|Ackermann et aL| 


2012 and the references therein). However, this is still well 
below PeV energies where CSF CRs are expected to be dom¬ 
inant and the Cherenkov Telescope Array | Actis et al. 
would be needed to constrain the gamma-ray spectra in PeV 
energy regime. More quantitative models of the CSF con¬ 
tribution to the diffuse 7 -ray and neutrino backgrounds will 
require better statistics of CSF SNe in starburst regions, ac¬ 
curate models of CR escape from the sources, as well as re¬ 
alistic models of CR propagation in starburst regions. 


6 CONCLUSIONS 

We have presented a colliding shock flow model for CR pro¬ 
duction in compact stellar clusters that efficiently produces 
hard CR spectra and neutrinos. Acceleration in colliding plas¬ 
mas is a very efficient version of Fermi acceleration given the 
strong confinement of CRs in the converging flows. The mech¬ 
anism is strongly nonlinear and time-dependent. We simulat¬ 
ed the acceleration process in a simplified geometry. Protons 
escape the accelerator with varying hardness and maximum 
energy as the SN shock approaches the wind. Furthermore, 
there is uncertainty in details of the mass distribution in the 
complicated outer ISM region (e.g., dense clouds) which will 
influence the results. Nevertheless, we believe our simulations 
include enough essential physics to estimate the neutrino and 
7 -ray emission from the galactic cluster Westerlund 1 and we 
show it is a likely source for IceCube events detected from 
the inner galaxy. Our 7 -ray predictions are consistent with 
H.E.S.S observations as well. 

While the relatively large angular uncertainty in the ar¬ 
rival directions of PeV neutrinos precludes an exact iden¬ 
tification, we believe some PeV IceCube events may result 
from ^ 10 PeV CR protons accelerated in Wdl (see Fig. [^. 
This cluster is a good candidate because it is one of the most 
massive clusters in the local group of galaxies and has an ob¬ 
served 10 ^ yr old magnetar, allowing enough time to spread 
PeV CRs over a few hundred parsecs scale. 

Future work that is critical for determining if a galac¬ 
tic supercluster can explain the apparent clump of 4-5 Ice- 
Cube events includes developing a more accurate model of 
multi-PeV CR diffusion on kpc scales. This requires a careful 
treatment of the matter distribution within a few degrees of 
Wdl and will result in a more accurate determination of the 
neutrino flux, as well as radio to 7 -ray emission, from CR 
interactions. 
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